% This function reads the ncfiles extracted from MITgcm model runs.
% This only reads the 109 level runs. To read 153 and 264 level runs
% read_ncfiles_km_runs.m should be used.

% Ritabrata Thakur, 18 Jan 2022.

function [data,lat,lon,ix,jx,ncpath,XC,YC,RC,time_start,time_end,time_array,z_ip] = read_ncfiles(vert_levs,MMP_number,run_no) 

basepath = <location of the ncfiles>;

%% section for 109 level
if (vert_levs == 109)
if (MMP_number == 1)
 if (run_no == 19)   
    ncpath = [basepath 'mmp_location/' 'lev_109_kpp_all_on_run19_mmp_loc' '.nc'];
 end
 if (run_no == 20)
    ncpath = [basepath 'mmp_location/' 'lev_109_kpp_back_visc_off_run20_mmp_loc' '.nc'];
 end
 lat = 25;
 lon = 195;
 ix = 90:91; %x-location of the model output
 jx = 87:88; %y-location of the model output
 % reading variables
data.temp = ncread(ncpath, 'T'); 
data.salt = ncread(ncpath, 'S');
data.U = ncread(ncpath, 'U');   
data.V = ncread(ncpath, 'V'); 
data.W = ncread(ncpath, 'W');
% you have to put here the path to the grids
gridx_name = 
gridy_name = 
gridz_name =  
XC = rdmds(gridx_name);
YC = rdmds(gridy_name);
RC = rdmds(gridz_name);
time_array = linspace(1, length(data.U(1,1,1,:)), length(data.U(1,1,1,:)));
time_start = time_array(1);
time_end   = time_array(end);
z_ip = RC;
end
% for MMP3
if (MMP_number == 3)    
 if (run_no == 19)   
    ncpath = [basepath 'mmp3_location/' 'lev_109_kpp_all_on_run19_mmp3_loc' '.nc'];
 end
 if (run_no == 20)
  ncpath = [basepath 'mmp3_location/' 'lev_109_kpp_back_visc_off_run20_mmp3_loc' '.nc'];
 end
 lat = 29;
 lon = 196.5;
 ix = 168:169; %x-location of the model output 
 jx = 288:289; %y-location of the model output
% reading variables
data.temp = ncread(ncpath, 'T'); 
data.salt = ncread(ncpath, 'S');
data.U = ncread(ncpath, 'U');   
data.V = ncread(ncpath, 'V'); 
data.W = ncread(ncpath, 'W');
% you have to put here the path to the grids
gridx_name = 
gridy_name = 
gridz_name = 
XC = rdmds(gridx_name);
YC = rdmds(gridy_name);
RC = rdmds(gridz_name);
time_array = linspace(1, length(data.U(1,1,1,:)), length(data.U(1,1,1,:)));
time_start = time_array(1);
time_end   = time_array(end);
z_ip = RC;
end
% for MP4
if (MMP_number == 4)
 if (run_no == 19)   
  ncpath = [basepath 'mmp4_location/' 'lev_109_kpp_all_on_run19_mmp4_loc' '.nc'];
 end
 if (run_no == 20)
  ncpath = [basepath 'mmp4_location/' 'lev_109_kpp_back_visc_off_run20_mmp4_loc' '.nc'];   
 end
 lat = 30;
 lon = 197;
 ix = 195:196; %x-location of the model output 
 jx = 360:361; %y-location of the model output
 % reading variables
data.temp = ncread(ncpath, 'T'); 
data.salt = ncread(ncpath, 'S');
data.U = ncread(ncpath, 'U');   
data.V = ncread(ncpath, 'V'); 
data.W = ncread(ncpath, 'W');
% you have to put here the path to the grids
gridx_name = 
gridy_name = 
gridz_name =
XC = rdmds(gridx_name);
YC = rdmds(gridy_name);
RC = rdmds(gridz_name);
time_array = linspace(1, length(data.U(1,1,1,:)), length(data.U(1,1,1,:)));
time_start = time_array(1);
time_end   = time_array(end);
z_ip = RC;
end
end

%% section for 153 level
if (vert_levs == 153)
if (MMP_number == 1)
 if (run_no == 34)   
    ncpath = [basepath 'mmp_location/' 'lev_153_kpp_all_on_run34_mmp1_loc' '.nc'];
 end
 if (run_no == 35)
    ncpath = [basepath 'mmp_location/' 'lev_153_back_off_run35_mmp1_loc' '.nc'];
 end
 lat = 25;
 lon = 195;
 ix = 90:91; %x-location of the model output 
 jx = 87:88; %y-location of the model output
  % reading variables
data.temp = ncread(ncpath, 'T'); 
data.salt = ncread(ncpath, 'S');
data.U = ncread(ncpath, 'U');   
data.V = ncread(ncpath, 'V'); 
data.W = ncread(ncpath, 'W');
% you have to put here the path to the grids
gridx_name = 
gridy_name = 
gridz_name = 
XC = rdmds(gridx_name);
YC = rdmds(gridy_name);
RC = rdmds(gridz_name);
time_array = linspace(1, length(data.U(1,1,1,:)), length(data.U(1,1,1,:)));
time_start = time_array(1);
time_end   = time_array(end);
z_ip = RC;
end
% for MMP3
if (MMP_number == 3)    
 if (run_no == 34)   
    ncpath = [basepath 'mmp3_location/' 'lev_153_kpp_all_on_run34_mmp3_loc' '.nc'];
 end
 if (run_no == 35)
  ncpath = [basepath 'mmp3_location/' 'lev_153_back_off_run35_mmp3_loc' '.nc'];
 end
 lat = 29;
 lon = 196.5;
 ix = 168:169; %x-location of the model output 
 jx = 288:289; %y-location of the model output
  % reading variables
data.temp = ncread(ncpath, 'T'); 
data.salt = ncread(ncpath, 'S');
data.U = ncread(ncpath, 'U');   
data.V = ncread(ncpath, 'V'); 
data.W = ncread(ncpath, 'W');
% you have to put here the path to the grids
gridx_name = 
gridy_name = 
gridz_name = 
XC = rdmds(gridx_name);
YC = rdmds(gridy_name);
RC = rdmds(gridz_name);
time_array = linspace(1, length(data.U(1,1,1,:)), length(data.U(1,1,1,:)));
time_start = time_array(1);
time_end   = time_array(end);
z_ip = RC;
end
% for MP4
if (MMP_number == 4)
 if (run_no == 34)   
  ncpath = [basepath 'mmp4_location/' 'lev_153_kpp_all_on_run34_mmp4_loc' '.nc'];
 end
 if (run_no == 35)
  ncpath = [basepath 'mmp4_location/' 'lev_153_back_off_run35_mmp4_loc' '.nc'];   
 end
 lat = 30;
 lon = 197;
 ix = 195:196; %x-location of the model output 
 jx = 360:361; %y-location of the model output
  % reading variables
data.temp = ncread(ncpath, 'T'); 
data.salt = ncread(ncpath, 'S');
data.U = ncread(ncpath, 'U');   
data.V = ncread(ncpath, 'V'); 
data.W = ncread(ncpath, 'W');
% you have to put here the path to the grids
gridx_name = 
gridy_name = 
gridz_name = 
XC = rdmds(gridx_name);
YC = rdmds(gridy_name);
RC = rdmds(gridz_name);
time_array = linspace(1, length(data.U(1,1,1,:)), length(data.U(1,1,1,:)));
time_start = time_array(1);
time_end   = time_array(end);
z_ip = RC;
end
end

%% section for 264 level
if (vert_levs == 264)
if (MMP_number == 1)
 if (run_no == 25)   
    ncpath = [basepath 'mmp_location/' 'lev_264_kpp_all_on_run25_mmp_loc' '.nc'];
 end
 if (run_no == 26)
    ncpath = [basepath 'mmp_location/' 'lev_264_kpp_back_visc_off_run26_mmp1_loc' '.nc'];
 end
 lat = 25;
 lon = 195;
 ix = 90:91; %x-location of the model output 
 jx = 87:88; %y-location of the model output
  % reading variables
data.temp = ncread(ncpath, 'T'); 
data.salt = ncread(ncpath, 'S');
data.U = ncread(ncpath, 'U');   
data.V = ncread(ncpath, 'V'); 
data.W = ncread(ncpath, 'W');
% you have to put here the path to the grids
gridx_name = 
gridy_name = 
gridz_name = 
XC = rdmds(gridx_name);
YC = rdmds(gridy_name);
RC = rdmds(gridz_name);
time_array = linspace(1, length(data.U(1,1,1,:)), length(data.U(1,1,1,:)));
time_start = time_array(1);
time_end   = time_array(end);
z_ip = RC;
end
% for MMP3
if (MMP_number == 3)    
 if (run_no == 25)   
    ncpath = [basepath 'mmp3_location/' 'lev_264_kpp_all_on_run25_mmp3_loc' '.nc'];
 end
 if (run_no == 26)
  ncpath = [basepath 'mmp3_location/' 'lev_264_kpp_back_visc_off_run26_mmp3_loc' '.nc'];
 end
 lat = 29;
 lon = 196.5;
 ix = 168:169; %x-location of the model output 
 jx = 288:289; %y-location of the model output
  % reading variables
data.temp = ncread(ncpath, 'T'); 
data.salt = ncread(ncpath, 'S');
data.U = ncread(ncpath, 'U');   
data.V = ncread(ncpath, 'V'); 
data.W = ncread(ncpath, 'W');
% you have to put here the path to the grids
gridx_name = 
gridy_name = 
gridz_name = 
XC = rdmds(gridx_name);
YC = rdmds(gridy_name);
RC = rdmds(gridz_name);
time_array = linspace(1, length(data.U(1,1,1,:)), length(data.U(1,1,1,:)));
time_start = time_array(1);
time_end   = time_array(end);
z_ip = RC;
end
% for MP4
if (MMP_number == 4)
 if (run_no == 25)   
  ncpath = [basepath 'mmp4_location/' 'lev_264_kpp_all_on_run25_mmp4_loc' '.nc'];
 end
 if (run_no == 26)
  ncpath = [basepath 'mmp4_location/' 'lev_264_kpp_back_visc_off_run26_mmp4_loc' '.nc'];   
 end
 lat = 30;
 lon = 197;
 ix = 195:196; %x-location of the model output 
 jx = 360:361; %y-location of the model output
  % reading variables
data.temp = ncread(ncpath, 'T'); 
data.salt = ncread(ncpath, 'S');
data.U = ncread(ncpath, 'U');   
data.V = ncread(ncpath, 'V'); 
data.W = ncread(ncpath, 'W');
% you have to put here the path to the grids
gridx_name = 
gridy_name = 
gridz_name = 
XC = rdmds(gridx_name);
YC = rdmds(gridy_name);
RC = rdmds(gridz_name);
time_array = linspace(1, length(data.U(1,1,1,:)), length(data.U(1,1,1,:)));
time_start = time_array(1);
time_end   = time_array(end);
z_ip = RC;
end
end
